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Development  of  a  3-D 
Tree  Thermal  Response  Model  for  Energy  Budget 

and  Scene  Simulation  Studies 


1  INTRODUCTION 

1.1  Background  and  Purpose  of  Research 

The  Balanced  Technology  Initiative  on  Smart  Weapons  Operability  Enhaixce- 
ment  (BTI/SWOE)  has  as  a  goal  the  ability  to  model  and  study  the  radiant  field 
from  complex  nahiral  backgrounds.  Remote  sensing  and  modeling  of  vegetated 
surfaces  is  a  complicated  endeavor  due  to  the  extreme  variability  that  can  be  en¬ 
countered  in  the  species  that  are  present.  The  sizes  of  vegetative  objects  are  highly 
variable  as  well  as  the  orientations  and  the  number  of  elements  and  components 
that  make  up  the  individual  tree  or  plant. 

The  purpose  of  this  study  has  been  to  develop  a  three  dimensional  thermal 
model  of  individual  trees  that  can  be  used  to  understand  the  thermal  properties 
of  trees.  The  model,  TREETHERM,  is  one  component  in  the  BTI/SWOE  Interim 
Thermal  Model  package,^  TREETHERM  is  designed  to  be  used  as  a  component  in 


Hummel,  J.R.,  Longrin,  D.R.,  Paul,  N.L.,  and  Jones,  J.R.  (1991)  “Development  of  the  Smart 
Weapons  Operability  Enhancement  Interim  Thermal  Model,”  Phillips  Laboratory,  Hanscom  AFB, 
Lexington,  Massachusetts,  PL-TR-91-2073,  1991. 
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remote  sensing  simulations  and  as  a  standalone  model  to  study  the  thermal  response 
characteristics  of  trees  to  changes  in  environmental  conditions. 

1.2  Organization  of  Report 

Section  2  describes  the  heat  conduction  model  that  was  developed  for  the  woody 
material  in  a  tree.  Section  3  describes  how  leaves  were  treated  in  the  model.  Sec¬ 
tion  4  presents  results  from  validation  efforts  made  with  the  model.  Finally,  Section 
5  presents  a  summary  of  our  efforts  and  recommendations  for  future  research.  A 
companion  report  provides  specific  details  on  the  use  of  TREETHERM  and  the 
data  required  to  run  the  code.'^ 


2  HEAT  CONDUCTION  MODEL  FOR  WOODY  MATERIAL 

2.1  General  Approach 

The  3-D  tree  model,  TREETHERM,  is  designed  to  be  used  in  scene  simulation 
studies  and  as  a  standalone  research  tool.  The  model  can  work  with  data  for 
specific  trees  or  be  used  to  model  generic  trees. 

In  scene  simulation  studies,  it  is  assumed  that  “standard”  trees  will  be  used 
for  given  tree  types.  Within  a  scene,  the  standard  tree  can  be  given  different 
orientations.  In  future  versions  of  the  model  it  is  envisioned  that  the  standard  tree 
could  be  scaled  to  different  sizes. 

2.2  Describing  the  Tree  Geometry 

The  modeling  philosophy  taken  in  developing  the  3-D  tree  model  has  been  to 
adapt  a  formalism  that  would  permit  the  user  to  model  “generic”  trees  or  specific 
trees.  To  achieve  the  latter  goal,  our  modeling  approach  has  been  to  base  our  3-D 
geometric  representation  of  trees  on  a  data  collection  approach  developed  by  the 
U.S.  Army  Corps  of  Engineers  Waterways  Experiment  Station  (WES)."" 

The  WES  database  describes  trees  by  giving  the  nodal  information  on  the 
locations  of  all  tree  components  and  the  starting,  ending,  and  average  diameters 
of  the  components.  The  nodal  information  is  the  x,  ij,  z  positions  relative  to  a 
reference  point  where  a  tree  element  either  changes  direction  or  connects  with 
another  element.  In  the  WES  data  base,  x  is  north,  y  is  west,  and  2  is  up.  A 

^  Jones,  J.R.  (1991)  “User’s  Guide  for  TREETHERM:  A  3-D  Thermal  Model  for  Single  Trees,” 
Phillips  Laboratory,  Hanscom  AFB,  Massachusetts,  PL-  rR-91-21()9,  31  March  1991. 

^  West,  H.W.,  and  Allen,  H.H.  (1971)  “A  Technique  for  Quantifying  Forest  Stands  for  Management 
Evaluations”,  US  Army  Engineer  Waterways  Experiment  Station,  Vicksburg,  MS,  Technical 
Report  M-71-9,  December  1971. 
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sampL  of  the  data  is  shown  in  Table  1.  A  prototype  user  interface  system  has 
also  been  developed  for  our  3-D  tree  thermal  model.  One  option  in  the  system 
is  to  access  the  database  and  display  a  graphical  representation  of  the  tree  being 
modeled.  Figure  1  shows  a  divSplay  from  a  prototype  user  interface  system  that 
was  developed  for  TREETHERM. 


HICKORY  TREE 
Texarkana,  Texas 


No.  of  Elenents:  326 
UaxlnuM  Height:  996  cb 
Trunk; 

#  elenents:  36 

Avg.  dianeter:  7.4  cm 
Branches : 

#  elenents:  290 

Avg.  dianeter:  1.3  cn 


*Fron  U.S.  Arnr  Engineer 
Hateruays  Experinent 
Station  Database 


Figure  1.  Graphical  Representation  of  a  Tree  as  Displayed  by  the  Prototype 
TREETHERM  User  Interface  System 


2.3  Developing  the  Computational  Framework 

Once  the  data  for  a  tree  have  been  assembled,  the  tree  components  are  modeled 
via  a  solid  modeling  technique  in  which  primitive  geometric  shapes  are  used  for 
each  component,  as  shown  in  Figure  2.  In  our  model,  truncated  cylinders  are  the 
shapes  being  used  to  represent  the  tree  components.  The  components  are  then 
divided  into  standardized  parts,  such  as  trunk,  limbs,  branches,  twigs,  etc  Once 
the  parts  are  identified,  the  material  distribution  within  each  part  is  established,  as 
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Table  1 .  Sample  of  the  Tree  Database  Collected  and  Assembled  by  the  Waterways 
Experiment  Station,  (The  number  trunk  and  branch  elements  was  determined  from 
the  data  and  is  not  a  part  of  the  WES  data.) 


Node 

# 

X 

Coord 

(cm) 

SOURCE 

Y  Z 

Coord  Coord 
(cm)  (cm) 

Diam 

(cm) 

TERMINUS 

X  Y  Z 

Node  Coord  Coord  Coord 
#  (cm)  (cm)  (cm) 

Diam 

(cm) 

Avg. 

Stem 

Diam 

(cm) 

1 

1379.0 

284.0 

13.0 

19.7 

2 

1376.0 

280.0 

81.0 

16.5 

18.1 

3 

1376.0 

280.0 

81.0 

16.5 

4 

1380.0 

287.0 

156.0 

12.0 

14,2 

e 

1380.0 

287.0 

156.0 

12.0 

6 

1375.0 

287.0 

231.0 

11.3 

11.6 

1 

1375.0 

287.0 

231.0 

1.4 

8 

1363.0 

311.0 

246.0 

0.7 

1.1 

9 

1375.0 

287.0 

231.0 

11.3 

10 

1375.0 

293.0 

299.0 

9.9 

10.6 

11 

1374.0 

293.0 

299.0 

2.5 

12 

1378.0 

306.0 

297.0 

1.7 

2.1 

13 

1378.0 

306.0 

297,0 

0.6 

14 

1377.0 

305.0 

300.0 

0.4 

0.5 

15 

1378.0 

306.0 

297.0 

0.6 

16 

1376.0 

305.0 

295.0 

0.4 

0.5 

17 

1378.0 

306.0 

297.0 

0.9 

18 

1375.0 

302.0 

297.0 

0.7 

0.8 

19 

1375.0 

293.0 

299.0 

1.7 

20 

1 374.0 

293.0 

317.0 

9.2 

5.4 

21 

1374.0 

293.0 

317.0 

2.8 

22 

1388.0 

307.0 

319.0 

2.5 

2.7 

23 

1388.0 

307.0 

319.0 

2.1 

24 

1336.0 

287.0 

321.0 

1.7 

1.9 

25 

1336.0 

287.0 

321.0 

1.7 

26 

1319.0 

283.0 

328.0 

1.5 

1.6 

27 

1319.0 

283.0 

328.0 

0.8 

28 

1315.0 

289.0 

330.0 

0.7 

0.7 

29 

1319.0 

283.0 

328.0 

1.5 

30 

1316.0 

283.0 

329.0 

1.5 

1.5 

31 

1316.0 

283.0 

329.0 

1.2 

32 

1300.0 

290.0 

338.0 

0.8 

1.0 

33 

1316.0 

283.0 

329.0 

1.5 

34 

1289.0 

294.0 

339.0 

0.5 

1.0 

35 

1370.0 

278.0 

316.0 

2.3 

36 

1349.0 

302.0 

317.0 

2.1 

2.2 

37 

1349.0 

302.0 

317.0 

2.1 

38 

1356.0 

337.0 

323.0 

1.8 

2.0 

39 

1356.0 

337.0 

323.0 

1.3 

40 

1345.0 

343.0 

325.0 

0.6 

0.9 

41 

1356.0 

337.0 

323.0 

1.4 

42 

1360.0 

361.0 

325.0 

0.9 

1.1 

43 

1374.0 

293.0 

317.0 

10.3 

44 

1378.0 

297.0 

332.0 

9.7 

10.0 

45 

1378.0 

297.0 

332.0 

2.6 

46 

1365.0 

260.0 

327.0 

2.4 

2.5 

47 

1365.0 

260.0 

327.0 

2.0 

48 

1367.0 

255.0 

329.0 

1.5 

1.8 

49 

1367.0 

255.0 

329.0 

0.8 

50 

1371.0 

260.0 

333.0 

0.6 

0.7 

shown  as  Step  1  in  Figure  3.  In  this  version  of  the  tree  model,  the  internal  material 
distribution  is  assumed  to  be  isotropic.  The  radii  of  the  cylinders  denoting  the 
different  woody  materials  are  user  inputs^  and  can  vary  from  tree  to  tree  and 
between  tree  components. 

Different  material  properties  can  be  assigned  to  the  different  tree  components. 
The  material  properties  can  vary  from  segment  to  segment.  This  feature  would 
allow  one  to  model  the  presence  of  dead  material  on  a  living  tree.  The  material 
properties  can  also  be  treated  as  temperature  dependent.  Within  each  modeled 
cylindrical  segment  the  material  distribution  is  assumed  to  be  homogeneous. 

This  version  does  not  treat  the  flow  of  moisture  within  the  tree.  However,  this 
can  be  examined,  at  least  to  a  first  approximation,  by  the  model’s  ability  to  vary 
material  properties  by  position. 
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Real  Tree 


Solid  Geometry 
Representation 


Figure  2.  Representing  Trees  as  a  Collection  of  Primitive  Solid  Geometric  Shapes 


Step  1  Step  2 

Establish  Internal  Divide  Up  Into 

Material  Distribution  Computational  Elements 


Step 
Formulate 
Energy 


Surface 

dget 


Bii 


Ollluis 

Solar 

FIuk 


Sanalbla 
^  Entrgy 
Flux 


IncORilma 
IR  Flux 


Latin  I 
Enargy 
Flux 


Step  4 

Perform  Lumped  Mass 
Calculation  With  Electrical 
Analog 


Figure  3.  Steps  Used  to  Perform  3-D  Thermal  Calculations  for  Trees 
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Each  tree  component  is  divided  into  cylindrical  subcomponents  and  each  sub¬ 
component  is  then  divided  into  computational  elements,  as  shown  as  Step  2  in 
Figure  3. 

2.4  Specifying  the  Surface  Energy  Budget 

The  surface  energy  budget,  as  shown  as  Step  3  in  Figure  3,  is  formulated  in 
terms  of  the  inconiing  solar  and  infrared,  reradiation  to  the  atmosphere,  convective 
and  latent  terms,  and  conduction  into  the  material.  Direct  and  diffuse  solar  fluxes 
are  attenuated  through  the  leafy  canopy  and  reflected  off  of  the  bark.  Self-shading 
from  the  various  tree  components  are  also  considered  via  the  use  of  a  ray-tracing 
routine.  The  tree  shading  model  uses  ray  casting  technology  initially  implemented 
in  SPARTA’S  Optical  Sensor  Simulation  (SENSORSIM)'*’^'®  for  rendering  images 
or  signatures  of  various  objects  as  observed  by  optical  sensor  systems.  The  con¬ 
vective  term  is  treated  in  terms  of  the  energy  transfer  resulting  from  air  flow  around 
cylindrical  surfaces. 


2.5  Temperature  Calculations 

The  temperatures  are  calculated  at  the  radial  mid-point  of  each  element  and 
applied  throughout  the  element.  The  heat  conduction  calculations  are  performed 
using  the  “lumped  mass”  approach  and  an  electrical  analog,  as  noted  as  Step  4  in 
Figure  3.  This  approach  is  similar  to  that  done  in  a  study  by  Derby  and  Gates.^  Each 
element  typically  interacts  with  up  to  six  adjacent  nodes. 

2.6  Thermal  Response 

The  thermal  response  of  the  model  will  be  in  three  spatial  dimensions  (x,  y,  z) 
predicted  from  the  general  heat  balance  equation  for  a  nonhomogeneous,  anisotropic 
material  is  given  by 

=  +  +  +  (i) 

where  p  is  the  density  of  the  woody  material  in  kg  m~^;  c  is  the  specific  heat  of  the 
material  in  J  kg~^  K~^;  T  is  the  temperature  in  K;  t  is  the  time  in  seconds;  kx,  ky, 

Guivens,  Jr.,  N.R.,  and  Henshaw,  P.D.  (1986)  “Laser  Radar  Sensor  Model,”  SPARTA,  Inc., 
Lexington,  MA,  SPARTA  LTR-86-08,27  June  1986. 

^  Guivens,  Jr.,  N.R.,  and  Henshaw,  P.D.  (1986)  “Laser  Radar  Sensor  Simulation,”  SPARTA,  Inc., 
Uxington,  MA,  SPARTA  LTR-87-011,  26  November  1986. 

®  Henshaw,  P.D.,  and  Guivens,  Jr.,  N.R.  (1987)  “Active/Passive  Multiwavelength  Imaging  Sensor 
Simulation,”  Laser  Radar  II,  R.  J.  Becherer  and  R.  C.  Harney,  Ed.,  Proc.  SPIE  783:84-90. 

^  Derby,  R.W.,  and  Gates,  D.M.,  (1966)  The  Temperature  of  Tree  Trunks  -  Calculated  and  Ob¬ 
served,  Amer.  J.  Bot.,  53:580-587. 
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and  kz  are  the  conductivities  in  the  a:,  y,  z  directions,  respectively,  in  W  m  ^  K 
and  Q{x,y,z,t)  is  the  net  summation  of  the  surface  energy  fluxes. 

In  this  initial  version  of  TREETHERM,  the  tree  elements  will  be  modeled 
by  the  lumped  mass  approach  for  cylindrical  segments.  Therefore,  the  x,y,z 
directions  above  will  correspond  to  the  radial,  circumferential,  and  longitudinal 
(axial)  directions,  respectively,  as  shown  in  Figure  4. 


Figure  4.  Schematic  Representation  of  the  Coordinate  System  Used  in  the  Heat 
Conduction  Calculations 


The  Crank-Nicholson  finite  difference  method  is  used  to  approximate  the  partial 
derivatives  in  the  energy  balance  equation.  This  approach  gives  more  accurate 
results  than  the  fully  implicit  or  explicit  methods.  In  general,  the  solution  scheme 
is  unconditionally  stable,  but  can  be  subject  to  oscillations.® 

The  coefficients  obtained  from  the  Crank-Nicholson  procedure  are  used  to  solve 
for  the  element  temperatures  using  a  modified  Seidel  iterative  scheme.  The  nu¬ 
merical  implementation  of  the  heat  balance  is  taken  from  Duncaner  al,^ 

The  finite  difference  form  of  Eqn.  1  is 


(2) 

where  i  is  the  number  of  the  element  being  solved  for,  j  is  the  number  of  the 


(t'-t 

0  =  (c.,- 

j^i  \  L 

®  Carnahan,  B.  Luther,  H.A.,  and  Wilkes  (1969)  Applied  Numerical  Methods.  Wiley  &  Sons,  New 
York,  p.451. 

®  Duncan,  T.C.,  Farr,  J.L.,  Wassel,  T.,  and  Curtis,  R.J.,  Satellite  Laser  Vulnerability  Model,  Ther¬ 
mal  Model  User’s  Guide,  Air  Force  Weapons  Laboratory,  (Software  Documentation). 
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element  in  thermal  contact  with  element  i,  n  is  the  total  number  of  elements  in 
thermal  contact  with  z,  A <  is  the  time  step  in  seconds,  V  is  the  volume  of  element 
j  in  m^,  Ti  and  Tj  are  the  previous  temperatures  in  K  for  elements  i  and  j,  T-  and 
Tj  are  the  temperatures  being  solved  for,  Cij  is  the  heat  conduction  parameter  in 

W  K“^,  and  is  the  total  surface  flux  at  element  i.  Thermal  contact  means  that 
energy  may  be  transferred  between  elements  by  conduction. 

The  heat  conduction  parameter  is  determined  from  an  electrical  analog,  the 
lumped  mass  approach,  and  given  in  terms  of  a  thermal  resistance,  Rij, 

=  7^  (3) 

where  the  thermal  resistance  is  given  by 

Rh  =  ^  (4) 

In  the  above,  the  thermal  resistance,  in  K  W“^,  is  evaluated  (see  Figure  5)  from 
the  center  or  radial  midpoint  of  element  i  to  the  surface  of  i  in  contact  with  the 
adjacent  element  j.  This  distance  is  Lij  and  is  given  in  meters,  ki  is  the  thermal 
conductivity  in  W  m~^  of  element  i  in  the  element  j  direction,  and  Aij  is 

the  average  cross-sectional  area  in  m^  between  the  calculation  location  of  i  to  the 
surface  of  i  in  contact  with  j.  The  lumped  mass  electrical  analog  has  the  advantage 
of  handling  arbitrary  shapes  and  element  connections. 


Figure  5.  Geometry  and  Properties  Used  in  Determining  the  Thermal  Resistance 
Between  Elements  i  and  j 

It  should  be  noted  that  the  current  version  of  TREETHERM  can  handle  tem¬ 
perature  dependent  properties  for  the  specific  heat,  conductivity,  emissivity,  and 
absorptivity.  This  version  does  not  attempt  a  mass  balance  due  to  moisture  trans¬ 
port  within  the  tree.  Also,  water  content  is  not  directly  considered.  Its  effect  can 
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be  approximated  by  changing  the  material  properties  used,  however.  Including  the 
effects  of  moisture  transport  is  planned  for  future  versions  of  the  model. 

The  implementation  of  the  energy  balance  (Eqn.  2)  has  the  advantage  of  al¬ 
lowing  each  element  being  modeled  to  have  different  boundary  fluxes  applied  to 
it.  For  the  current  version  of  TREETHERM,  these  boundary  fluxes  are  permitted 
only  at  surface  elements.  The  surface  conditions  currently  modeled  are  the  solar 
energy,  surface  reradiation  to  environment,  ant  surface  convection  due  to  wind. 

2.7  Solution  Scheme 

TREETHERM  uses  an  adaptive  time  step  in  its  solution  scheme.  The  time 
step  is  controlled  by  convergence  criteria  in  the  solution  scheme  and  r  ma’^^imum 
allowable  temperature  change  for  each  element  (1.5  K)  over  a  calculation  time 
interval. 

If  the  solution  for  a  time  step  does  not  converge  within  the  iteration  limit  (200 
steps),  the  model  retreats  to  the  start  of  the  time  interval,  resets  time  dependent 
parameters,  reduces  the  time  step,  and  repeats  the  calculation.  If  the  solution 
converges,  the  magnitude  of  each  element’s  temperature  change  is  compared  to 
the  maximum  allowed  temperature  change.  If  any  element’s  temperature  change 
over  the  time  step  exceeds  the  allowable,  the  calculation  will  be  repeated  with  the 
time  step  reduced.  This  time  step  reduction  for  the  repeat  calculation  will  take  into 
account  the  amount  by  which  the  temperature  change  criteria  was  exceeded.  If 
the  calculation  for  a  time  step  is  repeated  over  a  certain  number  of  intervals  (1(X) 
steps),  a  warning  is  issued  about  subsequent  results  and  the  calculation  allowed  to 
advance  in  time. 

If  the  solution  converges  and  the  element  temperature  change  criteria  are  met, 
the  model  updates  the  applicable  parameters  and  advances  in  time  by  the  current 
time  step.  The  scheme  will  try  to  increase  the  time  step  for  the  next  calculation 
based  on  the  ratio  of  the  maximum  element  temperature  change  to  the  allowable 
element  temperature  change  for  the  previous  time  step.  To  avoid  excess  repeat 
calculations,  the  maximum  increase  in  a  time  step  over  a  previous  time  step  is 
limited  to  10  per  cent.  A  time  step  is  also  subject  to  a  maximum  of  300  seconds 
and  a  minimum  of  1  second.  The  minimum  allowed  time  step  can  be  circumvented 
when  the  model  is  trying  to  hit  a  meteorological  data  time  interval. 

The  temperature  increase  criteria  and  the  maximum  and  minimum  allowed 
time  step  values  are  based  on  a  variety  of  factors.  These  factors  include  the 
expected  range  of  tree  material  properties  and  environmental  boundary  conditions, 
convergence  in  the  solution  scheme,  and  expected  meteorological  data  interval 
times,  as  well  as  exercising  the  code  for  various  combinations  of  temperature 
increase  and  time  step  criteria. 
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The  solution  iteration  and  repeat  calculation  limits  are  conservative  for  the  ex¬ 
pected  combination  of  tree  material  properties  and  environmental  boundary  condi¬ 
tions.  These  limits  are  applicable  to  the  general  purpose  three  dimensional  thermal 
code  from  which  TREETHERM  is  derived.  This  approach  is  also  used  in  the 
SWOETHRM  portion  of  the  BTI/SWOE  Interim  Thermal  Model  (ITM).^ 

2.7.1  Solar  Heat  Flux 

The  solar  heat  flux  is  a  time  and  spatial  varying  energy  source.  It  can  be 
obtained  from  data  or  estimated  from  various  radiative  transfer  parameterizations. 
The  approach  used  in  TREETHERM  is  the  same  as  that  used  in  the  Interim  Thermal 
Model.l 

The  solar  flux  incident  on  the  tree  elements  are  separated  into  direct  and  diffuse 
components.  The  diffuse  part  is  treated  as  radially  symmetric  to  a  given  tree 
segment.  The  total  solar  heat  flux  in  W  m~^  absorbed  by  a  surface  element, 
QHFi^  is  given  as 

QHFi  =  ai[Scir{x,ij,z,  t) A,  cos Sdfix,  y,  z,t) A, ^f]  (5) 

where  Sy^y.{x,y,  z,t)  and  Syf{x,  y,  z,t)  are  the  direct  and  diffuse  solar  energy, 
respectively,  on  element  i;  ai  is  the  surface  absorptivity  of  element  i;  Ai  is  the 
area  of  element  i  receiving  the  direct  solar  energy;  Ai^f  is  the  area  of  element  i 
receiving  the  diffuse  solar  energy;  and  9^  is  the  angle  of  incidence  between  the 
surface  normal  for  element  i  and  the  solar  vector. 

Shading  due  to  leaves  and  brtmches  as  well  as  attenuation  by  the  leaves  can 
modify  the  solar  energy  incidence  on  a  surface  element.  These  effects  are  handled 
via  an  adaptive  ray  casting  scheme  that  is  discussed  in  Appendix  A.  The  ray 
casting  scheme  determines  the  illuminated  area  of  branches  projected  onto  the 
solar  vector  line  of  sight  with  an  attenuation  modification  if  leaves  are  present 
in  the  model.  Section  3  of  this  report  details  the  modeling  of  the  leaf  energy 
transmission.  Figure  6  presents  some  of  the  possible  outcomes  of  a  set  of  rays  cast 
through  a  tree  model. 

The  shading/attenuation  effects  will  change  as  a  function  of  time  of  day.  How¬ 
ever,  within  TREETHERM  they  are  implemented  at  the  midpoint  of  the  time 
interval  for  the  meteorological  data  containing  the  solar  fluxes  rather  than  at  every 
computational  timestep.  This  is  done  primarily  for  computational  reasons.  The 
raycasting  approach  is  very  computationally  intensive  and  results  in  a  significant 
increase  in  computational  times.  The  approach  used  here  is  believed  to  be  accept¬ 
able  considering  the  uncertainties  in  tree  and  leaf  material  properties,  environmental 
conditions,  leaf  and  orientation  variations  due  to  wind  speed  changes,  etc. 
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Figure  6.  Shading  and  Solar  Attenuation  due  to  Branches  and  l^ves 

2.7.2  Radiation  to  Environment 

Currently,  the  tree  model  has  two  implementations  of  surface  re-radiation.  First, 
the  surface  elements  can  radiate  to  a  background  temperature.  This  method  is  used 
if  the  downward  infrared  data  are  not  available  from  the  meteorological  data  and 
the  built-in  infrared  parameterization  is  disabled.  In  this  case,  the  expression  for 
the  reradiation  is 

Qre.  =  (ri  -  i;^)  (6) 

where  QjiEi  is  the  infrared  radiation  to  the  background  environment  from  element 
i,  is  the  emissivity  of  the  surface  element  i,  Ai  is  the  area  in  vc?  of  the  radiating 
surface,  SFi  is  the  shape  factor  (which  is  1  in  this  version),  a  is  the  Stefan- 
Boltzmann  constant,  and  Too  is  the  background  temperature  in  K  obtained  from 
the  meteorological  data.  The  second  method  of  determining  the  net  re-radiation  at 
a  surface  is 

Qrei  =  Oiij^AilR  -  e^AiSFioTf  (7) 

where  IR  is  the  downwelling  infrared  flux  from  the  atmosphere  in  W  m~^. 

The  infrared  flux  from  the  ground  is  not  modeled  here.  An  approximation 
to  this  is  done  by  letting  the  downward  atmospheric  infrared  flux  be  spherically 
symmetric.  Then,  adjustments  can  be  made  to  the  material  infrared  surface  ab¬ 
sorptivity.  That  is,  those  computational  elements  that  face  the  ground  can  have 
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an  infrared  surface  absorptivity  different  than  those  facing  the  sky.  In  future  ver¬ 
sions  of  TREETHERM,  the  upwelling  infrared  flux  from  a  thermal  model  of  the 
underlying  ground  will  be  added. 


2.7.3  Surface  Convection 

The  heat  transfer  between  the  surrounding  air  and  a  given  surface  element  of 
a  tree  segment  is  given  as 


QsCi  =  -  T,)  (8) 

where  QsCi  surface  convection  at  element  i,  hj  is  the  convective  coefficient 
at  element  ^  in  W  m“^  K“^,  and  Ai  is  the  surface  area  of  element  i  in  m^. 

The  coefficient  for  forced  convection,  which  occurs  when  the  wind  speed  is  not 
zero,  is  determined  from  Churchill’s  correlation  equations  for  cylinders  as^®’^^ 


If  ^  /l  / 

~D 


(9) 


where  Nuf)  is  the  mean  Nusselt  number,  kj  is  the  fluid  (air)  conductivity  in 

W  m“^  K~^,  and  D  is  the  cylinder  diameter  in  m.  The  mean  Nusselt  number  is 
determined  as  a  function  of  the  Reynolds  number,  Rej)\  Prandtl  number,  Pr;  or 
Peclet  number,  Pe.  For  example. 


7r^  =  0.3+ 

*  2 '  ^ 
l-i-(0.4/Pr)3 

=  0.3  -f 

[l+(0.4/Pr)i  ^ 


ReD<  4000  (10) 


+ 


20,000  <  <  400,000  (11) 


Nu£)  —  0.3  -f- 


1  1 
0.62Ren'^Pr^ 

l-|-(0.4/Pr)3 


General  Case  (12) 


Churchill,  S.A.,  and  Bernstein,  M.  (1977),  A  Correlating  Equation  for  Forced  Convection  from 
Gases  and  Liquid  to  a  Circular  Cylinder  in  Crossflow,  J.  Heat  Transfer,  99:300-306. 

Lienhard,  J.  M.,  (1981),  A  Heat  Transfer  Textbook.  Prentice-Hall,  Englewood  Clirr»,  NJ,  ISSl, 
pp  326-337. 
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Nud  = 


1 


Pe  <  0.2 


(13) 


0.8237-ln^PeJ 

The  expressions  for  the  Nusselt  number  given  above  are  based  on  empirical  obser¬ 
vations.  (The  reader  is  referred  to  any  standard  text  on  heat  transfer  for  a  discussion 
of  the  above  parameters.)  The  Reynolds  number  is  given  as 

Rcf)  =  uDjv  (14) 

where  u  is  the  free  stream  (wind)  velocity  in  m/s,  D  is  the  diameter  of  the  cylinder 
in  m,  and  v  is  the  fluid  (air)  kinematic  viscosity  in  m^  s~^.  The  Prandtl  number 
is  given  as 

Pr  =  vja  (15) 

where  ot  is  the  fluid  (air)  diffusivity  in  m^  s“^.  Finally,  the  Peclet  number  is  given 
as 


Pc  r-  Pr  .  Pe£).  (16) 

Eqn.  12  is  given  as  the  correlating  equation  over  the  entire  range  of  Rej), 
and  Eqns.  10  and  11  are  used  for  the  typical  applications  used  in  TREETHERM. 
Eqn.  13  is  recommended  when  the  Peclet  number  is  less  than  0.2. 

It  is  recognized  that  the  forced  convective  term  hi  will  vary  circumferentially. 
Giedt^^  has  provided  a  range  of  experimental  data  to  estimate  hi  as  a  function 
of  radial  location  from  the  stagnation  point.  These  data  were  used  in  a  study  of 
temperature  variation  in  tree  trunks^  and  are  similar  to  those  calculated  by  the 
model.  This  method  will  be  included  as  an  alternative  to  the  current  convective 
estimation  in  a  subsequent  version  of  the  tree  model. 

A  free  convective  term  has  not  yet  been  incorporated  into  the  convective  heat 
transfer  term.  However,  a  modification  factor  can  be  applied  to  the  Nusselt  number 
to  account  for  a  small  amount  of  assisting  or  opposing  free  convection. 


Giedt,  W.  H.,  (1949),  Investigation  of  Variation  of  Point  Unit  Heat-Transfer  Coefficient  Around 
a  Cylinder  Normal  to  an  Air  Stream,  Trans.  ASME,  71:375-381. 

Churchill,  S.  W.,  (1977),  A  Comprehensive  Correlating  Equation  for  Laminar,  Assisting,  Forced 
and  Free  Convection,  AIChE  Journal,  23:10-16. 
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3  TREATMENT  OF  LEAVES 

The  tree  geometry  database  used  to  describe  a  given  tree  can  contain  the  po¬ 
sitions  of  leaf  clusters  on  branches.  Currently,  we  are  only  considering  deciduous 
trees.  Future  versions  of  TREETHERM  will  be  expanded  to  include  coniferous 
trees. 

We  are  assuming  that  the  leaves  occur  in  cylindrical  clusters  around  the  branch 
segments,  as  shown  in  Figure  7.  We  are  also  assuming  that  the  leaves  are  uniform 
in  size  within  each  cluster,  the  clusters  have  a  medium  density  of  leaves,  and 
that  the  leaves  are  well  ventilated.  In  this  version  of  TREETHERM,  the  leaves 
are  treated  as  flat  plates  and  are  assumed  to  be  randomly  oriented  relative  to  tree 
elements  for  all  wind  speeds.  This  will  be  examined  in  later  versions  of  the  model, 
however. 


Figure  7.  Schematic  Representation  of  the  Leaf  Clusters 
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3.1  Leaf  Attenuation 

The  attenuation  of  radiation  by  the  leaf  clusters  is  being  treated  with  an  aerosol 
attenuation  analog.  Specifically,  the  leaf  attenuation  model  calculates  macroscopic 
attenuation  properties  by  summing  the  attenuation  effects  from  individual  leaves 
making  up  the  cluster.  The  macroscopic  attenuation  properties  account  for  gaps  in 
leaf  clusters  and  leaf  overlapping  as  described  below. 

The  general  approach  is  to  model  a  leaf  cluster  as  a  cylindrical  envelope  con¬ 
taining  leaves  and  air  where  the  parent  branch  segment  is  the  symmetry  axis.  The 
observed  geometry  of  individual  leaves  is  then  used  to  construct  macroscopic  val¬ 
ues  of  transmittance,  Tc,  absorptance,  Ac,  and  reflectance,  for  tho  loaf  cluster. 
Physically,  these  macroscopic  values  represent  the  fraction  of  energy  transmitted 
through,  absorbed  in,  and  reflected  by  the  leaf  cluster  as  a  whole.  Additionally, 
the  macroscopic  transmission  to  the  parent  branch,  is  also  calculated.  Any  ray 
encountering  a  leaf  cluster  will  interact  according  to  these  macroscopic  amounts. 
(More  precise  modeling  would  require  detailed  information  about  the  internal  stmc- 
ture  of  leaf  clusters  which  is  beyond  the  scope  of  this  effort.)  The  macroscopic 
parameters  of  the  leaf  clusters  have  the  properties 

Tc{0)lo  =  Radiation  Transmitted  Through  A  Cluster 

Ac{0)lo  =  Radiation  Absorbed  In  A  Clust<*r 

Rc{^)^o  =  Radiation  Reflected  By  A  Cluster 

T})c{0)lo  =  Radiation  Transmitted  To  A  Parent  Branch 

where  Iq  is  the  incident  radiation  and  0  is  the  angle  between  the  parent  branch 
segment  and  the  incident  ray.  Due  to  symmetry  arguments,  the  macroscopic  pa¬ 
rameters  have  the  property  that  Tc{0)  =  Tc(180  -  6),  for  example.  The  angular 
dependence  arises  because  the  potential  shading  area  of  a  cluster  depends  on  the 
direction  of  the  incident  radiation.  When  the  incident  ray  is  along  the  symmetry 
axis,  for  example,  leaves  cast  theii  shadow  into  the  area  defined  by  the  circular 
face  of  the  cylinder.  Similarly,  when  the  incident  ray  is  normal  to  the  symmetry 
axis,  leaf  shadows  are  in  a  rectangular  area  defined  by  the  diameter  and  length  of 
the  cylinder.  Formally,  the  potential  shading  area  of  a  cylindrical  cluster,  Xc  in 
cm^,  is 

Xc{9)  =  cos  0  +  2r  Ism  0  (17) 

where  r  and  I  are  the  radius  and  length  of  the  cluster,  respectively.  For  Tj^j.  how¬ 
ever,  it  is  necessary  to  describe  the  radiation  passing  through  the  plane  containing 
the  parent  branch.  Here,  the  potential  shading  area  is  the  projection  of  the  rectangle 
defined  by 
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Xi,j.{0)  =  2rlsmd.  (18) 

As  a  result  of  the  assumption  of  symmetry,  there  is  no  azimuthal  dependence  on 
either  Xc  or 

To  develop  macroscopic  leaf  parameters,  individual  aspen  leaves,  for  example, 
are  assumed  to  be  uniform  circular  plates  each  having  a  surface  area  of  30  cm“.  The 
leaf  number  density  is  taken  as  30  leaves  per  foot  of  leaf  cluster.  (In  TT^EETHERM, 
the  surface  area  of  individual  leaves  and  the  leaf  number  density  are  user  supplied 
inputs.)  Due  to  the  strong  influence  of  wind,  '  is  assumed  that  individual  leaves 
are  randomly  oriented  inside  the  cylindrical  cluster.  Speciflcally,  the  shadowed 
area  casted  by  an  individual  leaf,  Xi{(f>)  in  cm“,  is 

Xi{0)  =  30  cos  0  (19) 


where  0,  which  randomly  varies  between  0  and  90  degrees,  is  the  angle  between 
the  incident  radiation  and  tlie  normal  to  the  leaf  plate.  Partial  leaf  overlapping  is 
also  addressed  by  means  of  a  statistical  approach.  To  do  this,  individual  leaves  are 
randomly  placed  on  Xc{0)  such  that  the  total  area  shaded  by  leaves  is 


SciO)  = 


A1(0) 

oiXj{(p) 


e,  >  sr\e)/Xc(0) 

€i  <  sr'w/-Vc(0) 


(20) 


where  ii  is  the  number  of  leaves,  cj  and  are  random  numbers  between  0  and 
1,  and  5c“^(^)  /Xc{0)  is  the  cumulative  fractional  shading  up  to  the  i  -  1  leaf. 
Implicit  in  the  above  expression  is  that  leaves  can  be  placed  anywhere  within  the 
cylindrical  envelop.  Physically,  this  is  a  reasonable  assumption  given  that  1.)  leaf 
stems  vary  in  length  and  2.)  leaves  and  their  stems  are  often  bent  such  they  are 
close  to  the  parent  branch  during  windy  conditions. 

It  also  follows  that  the  area  covered  by  overlapping  leaves  is  given  by 


Oc{0)  =  E"=i(l  -  o,)Xj{4>)  C,  <  sr‘(9)/.v,(«). 


(21) 


Future  leaf  modeling  efforts  will  modify  Eqn  19  and  21  to  account  for  diffraction 
effects  around  leaf  edges  which  tend  to  blur  the  shadow  casted  by  an  individual 
leaf. 

The  above  expressions  for  Sc{0)  and  Oc{0)  have  been  evaluated  for  a  leaf 
cluster  4  feet  long  by  1  foot  wide.  Figure  8  shows  these  results  as  a  function 
of  incident  angle  (6)  where  Sc{0)  and  Oc{(f)  have  been  normalized  by 
The  model  suggests  that  about  60  percent  of  A',  is  shaded  when  incident  rays 
encounter  the  cluster  lengthwise  (40  <  0  <  90  deg)  which  reasonably  agrees  with 
visual  observations  made  of  an  aspen  used  in  our  validation  studies  (see  Section  4.) 
The  amount  of  shading  and  overlapping  increases  when  incident  angles  decrease 
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Incident  Angle  (deg) 


Figure  8.  Fraction  of  AV  Shaded  by  Leaves  as  a  Function  of  Incident  Angle  for 
a  4  Foot  Long  by  I  Foot  Wide  Leaf  Cluster.  The  fraction  shading  by  overlapping 
leaves  is  also  shown  in  the  Figure 


from  40°  to  20°.  Below  about  20°,  where  the  rays  are  directed  along  the  parent 
branch  axis,  Xc  is  completely  shadowed  by  leaves. 

The  above  expressions  for  leaf  shading  can  be  used  to  obtain  the  macroscopic 
transmission  coefficient  for  a  leaf  cluster 

^  ,n^  _  {^^r{O)-SAe)+ti(Sc(0)-(MO))+t‘fOcm 

(22) 

where  ti  is  the  transmission  for  an  individual  leaf.  (Angular  variations  in  t;  are  not 
considered  in  this  effort.)  Effectively,  the  individual  terms  in  the  above  expression 
represent  weighted  averages  of  the  areas  covered  by  no  leaves,  single  leaves  and 
overlapping  leaves.  The  overlapping  leaf  term  assumes  incident  rays  encounter 
two  leaves  before  exiting  which  may  not  be  valid  for  dense  leaf  clusters.  Also  the 
shading  effects  of  the  parent  branch  are  assumed  negligible.  Similar  expressions 
can  be  developed  for  the  macroscopic  absorption  and  reflection  coefficients 


i  in\  _  {<HSr{0)-k-aiiiOr{6)) 

~  XA^) 


(23) 


P  (  a\  _  (r^5r(^)-f 

-  XAO) 


[0) 

where  and  r/  are  the  absorption  and  reflection  for  an  individual  leaf. 


(24) 

The 
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macroscopic  parameters  have  the  property  Tc{9)  +  AdO)  +  Rc{9)  —  1. 

The  radiation  transmitted  to  a  parent  branch,  is  given  by 

rj.  _  {Xkr{e)-sud)-vti{sad)-o^rm+tfo^,[e)) 

lhr\9)  —  AV(^) 

where  is  given  by  Eqn  18,  5^,.  is  the  leaf  shading  in  the  plane  containing 
the  parent  branch,  and  is  the  area  covered  by  overlapping  leaves  in  the  plane 
containing  the  parent  branch.  The  leaf  shading.S^,.,  is  given  by 


= srii 


x\(o) 

o,-V|(<i) 


=.  >  Sl;\e) ! 

! Xy(e) 


and  Oi,,  is  given  by  a  form  similar  to  Eqn  21. 


(26) 


3.2  Calculating  Leaf  Temperatures 

3.2.1  Leaf  Energy  Budget 

The  leaf  energy  budget  model  is  based  on  a  formulation  of  Gates'^  for  a  single 
leaf.  The  energy  budget  is  given  as 

a,t,.(S  +  s)+ajfrs{S  +  3)  ^{Ra  +  Rg)  _  ^^^^4  ±  C  ±  LE  =  0  (27) 

where  is  the  leaf  absorpiance  to  direct  solar  radiation,  a^j  is  the  leaf  absorp- 
tance  to  diffuse  solar  radiation,  S  is  the  direct  solar  flux,  s  is  the  diffuse  solar  flux, 
/’5  is  the  surface  albedo,  at  is  the  leaf  infrared  absorptance,  Ba  is  the  downwelling 
infrared  radiation  from  the  atmosphere,  Rg  is  the  upwelling  infrared  radiation  from 
the  underlying  ground,  (f  is  the  leaf  infrared  emittance,  a  the  Stefan-Boltzmann 
constant,  T[  is  the  leaf  temperature,  C  is  the  convective  exchange  term,  L  is  the 
latent  heat  of  vaporization,  and  E  is  the  evapotranspiration  term. 

The  convective  exchange  term  is  given  as^"^ 

C  =  hdT(-Ta)  (28) 

where  he  is  the  convection  coefficient  and  Ta  is  the  air  temperature.  The  convection 
coefficient  is  dependent  upon  the  shape  orientation,  and  size  of  the  leaf;  the  state 
of  the  air  flow  (laminar  or  turbulent);  whether  the  convection  is  free  or  forced; 
and  upon  any  interactions  between  leaves.  For  the  purposes  of  this  study,  the 
leaves  will  be  treated  as  circular  flat  plates  that  have  minimal  interactions  with  one 
another  and  are  oriented  into  the  wind.  In  the  case  of  free  convection  in  still  air, 
h(  can  be  given  as^'^ 

'  Gates,  D  M.  (1964)  “Characteristics  of  Soil  and  Vegetated  Surfaces  to  Reflected  and  Emitted 
Radiation,”  Proceedings  IGARD  Symposium  on  Remote  Sensing  of  Environment,  573-600. 
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/!c  =  4.19[3^P 


(29) 


where  is  the  leaf  diameter.  For  forced  convection,  he  is  given  as^^ 

hc  =  ^M[^y  (30) 

where  V  is  the  wind  speed. 

The  evaporation  of  water  from  a  leaf  occurs  as  a  process  of  water  vapor  dif¬ 
fusion  from  the  saturated  inner  cells  through  the  leaf  stomata  and  across  the  leaf 
boundary  layer  to  the  free  air.  The  rate  of  water  vapor  diffusion,  E,  is  proportional 
to  the  gradient  of  water  vapor  density  between  the  leaf  and  the  free  air.  This  can 
be  expressed  as 


where  pg  is  the  saturation  water  vapor  density  at  the  leaf  temperature,  rh  is 
the  relative  humidity  of  the  ambient  air,  ps{Ta)  is  the  saturation  water  vapor 
density  of  the  air,  and  R  is  the  diffusion  resistance  of  the  leaf.  The  diffusion 
resistance  depends  upon  the  leaf  morphology  of  the  given  species.  Gates^^  gives 
a  characteristic  value  of  6.0  sec  cm”^  for  many  deciduous  species  in  full  sunlight 
with  the  stomates  wide  open.  In  contrast,  the  diffusion  resistance  for  conifers  is 
given  by  Gates  as  30.0  sec  cm~^  or  greater.  The  diffusion  resistance  should  be 
treated  as  temporally  dependent  but,  for  the  purposes  of  this  first  generation  model, 
will  be  treated  as  a  constant. 

3.2.2  Impact  of  Model  Assumptions 

Representative  leaf  temperatures  were  calculated  using  weather  conditions  at 
Hunter-Liggett,  California,  on  20  September  1989.  This  was  the  day  chosen  for 
the  FY  90  BTI/SWOE  demonstration.  The  day  was  a  clear  day.  Table  2  lists  the 
weather  data  for  the  day  and  Table  3  lists  the  default  parameters  used  in  the  leaf 
temperature  calculations.  The  surface  albedo  is  representative  of  bare  surfaces  or 
surfaces  with  light  vegetation,  The  leaf  solar  absorptances  are  representative  of 
a  variety  of  plants.^® 


Hummel,  J.R.,  and  Reck,  R.A.  (1979)  A  Ulobal  Surface  Albedo  Model,  J.  Appl.  Meteor., 
18:239-253. 

Gates,  D.M.,  Kegan,  H.J.,  Schleter,  J.C.,  and  Weidner,  V.R.  (1965)  Spectral  Properties  of  Plants, 
Appl.Optics,  4:11-20. 


Table  2.  Weather  Data  for  the  FY  90  SWOE  Hunter-Liggett  Demonstration.  The 
data  are  for  Julian  Day  263,  20  September  1989 


TIME  TEMP 
(LST)  (C) 

CLOUD  WIND  WIND 

RH  VISIBILITY  CLOUD  COVER  SPEED  DIRECTION 
(%)  (km)  OBS  (%)  (m/iec)  (deg) 

SURFACE 

PRESS 

(mb) 

0100 

10.6 

89 

25 

Clear 

0 

0.31 

128 

968.8 

0200 

10.5 

90 

25 

Clear 

0 

0.67 

76 

969.2 

0300 

9.4 

90 

25 

Clear 

0 

0.05 

Calm 

969.2 

0400 

9.5 

91 

25 

Clear 

0 

0.21 

91 

969.4 

0500 

8.7 

91 

25 

Clear 

0 

0.62 

56 

969.6 

0600 

8.2 

91 

16 

St 

100 

0.57 

65 

969.9 

0700 

10.8 

91 

2 

Fog 

100 

0.10 

Calm 

970.6 

0800 

13.0 

87 

2 

Fog 

100 

0.26 

246 

970.8 

0900 

15.2 

77 

16 

St 

50 

0.98 

17 

971.1 

1000 

19.1 

55 

25 

Clear 

0 

1.34 

93 

971.1 

1100 

21.3 

43 

25 

Clear 

0 

1.59 

67 

970.6 

1200 

23.7 

34 

25 

Clear 

0 

1.39 

31 

970.1 

1300 

25.2 

21 

25 

Clear 

0 

1.95 

89 

969.6 

1400 

26.7 

21 

25 

Clear 

0 

1.49 

351 

969.0 

1500 

27.6 

18 

25 

Clear 

0 

1.39 

59 

968.5 

1600 

27.4 

17 

25 

Clear 

0 

2.36 

57 

968.1 

1700 

26.1 

20 

25 

Clear 

0 

3.03 

49 

968.3 

1800 

24.3 

23 

25 

Clear 

0 

1.18 

50 

968.3 

1900 

19.6 

41 

25 

Clear 

0 

0.26 

220 

968.5 

2000 

18.5 

45 

25 

Clear 

0 

1.13 

137 

969.0 

2100 

16.3 

54 

25 

Clear 

0 

0.72 

162 

969.3 

2200 

15.2 

59 

25 

Clear 

0 

0.51 

105 

969.5 

2300 

12.7 

75 

25 

Clear 

0 

0.36 

36 

970.0 

2400 

12.3 

81 

25 

Clear 

0 

0.21 

Calm 

970.0 

20 


Table  3.  Model  Conditions  Assumed  for  Leaf  Temperature  Calculations 


PARAMETER 

VALUE 

Surface  Albedo 

0.15 

Leaf  Width 

5  cm 

Diffusion  Resistance 

6  cm/sec 

Leaf  Absorptance 
Direct  Solar 

0.6 

Diffuse  Solar 

0.6 

Leaf  IR  Absorptivity 

0.97 

Leaf  IR  Emissivity 

0.97 

3.2.2.1  Impact  of  Evapotranspiration 

Figure  9  displays  the  leaf  temperatures  with  and  without  evapotranspiration 
considered.  The  air  temperature  is  also  shown  for  reference.  The  Figure  shows 
that,  in  general,  leaf  temperatures  are  lower  with  evapotranspiration  considered 
than  with  it  neglected.  During  the  period  of  peak  solar  loading,  the  neglect  of 
evapotranspiration  can  lead  to  differences  in  the  calculated  leaf  temperature  of  as 
high  as  eight  degrees. 


Figure  9.  Leaf  Temperatures  Calculated  With  and  Without  Evapotranspiration  for 
Weather  Conditions  Representative  of  the  FY  90  SWOE  Hunter-Liggett  Demon¬ 
stration 
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Figure  10.  Leaf  Temperatures  Calculated  With  Free  and  Forced  Convection  As¬ 
sumed  for  Weather  Conditions  Representative  of  the  FY  90  SWOE  Hunter-Liggett 
Demonstration 

3.2.2.2  Impact  of  Convection 

Figure  10  displays  the  leaf  temperatures  with  free  and  forced  convection  con¬ 
sidered.  The  air  temperature  is  also  shown  for  reference.  The  free  convection 
results  were  obtained  by  setting  the  wind  speeds  to  zero.  With  free  convection 
assumed,  the  leaf  temperature  generally  follows  the  variations  in  solar  loading. 
With  forced  convection  assumed,  the  leaf  temperature  also  varies  in  response  to 
the  wind  speed. 

3.2.2.3  Impact  of  Leaf  Diameter 

Figure  11  displays  the  leaf  temperatures  with  three  different  leaf  diameters 
considered.  The  results  show  that  the  leaf  temperature  is  relatively  insensitive  to 
the  assumed  leaf  diameter. 

4  VALIDATION  EFFORTS 
4.1  Overview  of  STAMP 

To  aid  in  the  development  of  the  3-D  tree  model,  SPARTA  designed  and  con¬ 
ducted  a  measurement  program  that  was  held  near  Bangor,  Maine.  The  purpose  of 
this  field  test  program  was  to  provide  a  database  for  validation  efforts  and  to  aid 
in  identifying  important  phenomenology. 

A  bigtooth  aspen  on  land  owned  by  the  International  Paper  Company  was 
selected  for  study.  The  measurement  program  was  held  6  -  24  September  1990 
and  occurred  in  conjunction  with  a  NASA  program  on  Forest  Ecosystem  Dynamics. 
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Local  Time  (hours) 

Figure  1 1 .  Leaf  Temperatures  Calculated  With  Different  Leaf  Diameters  Assumed 
for  Weather  Conditions  Representative  of  the  FY  90  SWOE  Hunter-Liggett  Demon¬ 
stration 

The  aspen  tree,  nicknamed  Joyce  Kilmer  (after  the  poet),  was  instrumented 
with  thermocouples  and  thermistors  to  measure  the  temperature  at  various  locations 
within  the  tree.  Leaf  and  bark  samples  were  also  taken  and  optical  measurements 
made.^^  Supporting  weather  data  were  provided  by  the  University  of  Maine  and 
imagery  data  were  taken  by  personnel  from  the  Keweenaw  Research  Center  of  the 
Michigan  Technical  University. 

4.2  Discussion  of  Results 

In  this  section,  the  TREETHERM  model  will  be  used  to  predict  temperatures 
at  various  locations  in  a  tree  trunk.  Results  will  be  given  for  comparison  against 
the  Maine  tree,  “Joyce  Kilmer”,  and  a  tree  located  at  Hunter-Liggett,  California. 
This  was  the  site  being  modeled  for  the  demonstration  effort  in  the  second  year  of 
SWOE. 


Neu,  J.T.,  Dummer,  R.S.,  Beecroft,  M.,  McKenna,  and  Robertson,  D.C.  (1990)  “Surface  Optical 
Property  Measurements  on  bark  and  Leaf  Samples,”  Geophysics  Laboratory,  Hanscom  AFB, 
Massachusetts,  PL-TR-9 1-2009 
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4.2.1  Temperature  Predictions  for  “Joyce  Kilmer” 

A  cross-section  of  the  trunk  4.5”  in  diameter  at  the  one  foot  level  was 
modeled.  Little  to  no  shading  occurred  from  other  branches  or  leaves  at  this 
location.  Meteorological  data  from  8  and  9  September  1990  were  used  to  calculate 
the  solar  and  downward  infrared  boundary  conditions.  The  meteorological  data 
were  taken  at  a  site  immediately  adjacent  to  the  tree.  Weather  conditions  were 
generally  clear  on  those  days  with  daytime  maximum  temperatures  near  20  C  and 
nighttime  minimum  temperatures  between  1  -  5  C. 

Temperature  predictions  from  TREETHERM  have  been  compared  against  ther¬ 
mocouple  data  taken  at  locations  on  the  north  and  south  faces  at  about  1”  in  from 
the  bark  surface.  Comparisons  for  two  days  of  the  model  results  and  measured 
values  are  shown  in  Figures  12  and  13,  respectively.  Additionally,  Figures  14  and 
15  show  the  predicted  temperature  cross  section  distributions  at  9  AM  and  3  PM 
for  9  September  1990.  In  general,  the  agreement  between  the  model  and  data  are 
good  but  some  nodel  limitations  are  apparent.  For  example,  the  peak  tree  temper¬ 
ature  is  overpredicted  by  the  model  which  then  leads  to  a  consistent  overprediction 
during  the  rest  of  the  day.  This  overprediction  could  be  the  result  of  a  number 
of  factors  such  as  incorrect  material  properties,  the  inability  to  account  for  water 
vapor  condensing  on  the  bark  (a  process  that  occurred  during  the  early  morning), 
and  uncertainties  in  the  solar  and  infrared  fluxes. 

The  average  convective  term  implemented  does  not  predict  the  temperature 
plateaus  from  the  south  face  thermocouple  data  on  both  8  and  9  September  1990. 
On  8  September  1990,  the  wind  during  the  period  9  AM  to  6  PM  varied  from 
easterly  to  southerly.  For  the  same  period  on  9  September  1990,  the  wind  was 
mostly  southerly.  The  thermocouple  data  shows  this  by  the  leveling  off  of  the 
temperatures  at  the  south  face  thermocouple  location.  The  wind  effects  are  not  as 
pronounced  in  the  north  face  thermocouple  data  as  that  location  was  generally  in 
the  lee  of  the  wind.  Modeling  the  effect  of  the  wind  direction  would  improve  the 
modeling  of  the  forced  convective  term  for  exposed  tmnks  and  branches.  Refine¬ 
ment  of  the  forced  convective  term  further  than  accounting  for  the  wind  direction 
and  possibly  the  angle  of  incidence  is  not  contemplated.  Any  effects  due  to  the 
three  dimensional  disturbances  in  the  velocity  field  are  beyond  the  scope  of  this 
model. 

The  model  also  did  not  predict  the  abrupt  morning  temperature  rise  very  well 
on  8  September  1990  but  did  a  better  job  on  9  September  1990.  This  is  most 
likely  a  consequence  of  the  choice  of  material  properties  and  the  estimation  of  the 
incident  energy  on  the  model  surface. 

Comparison  of  the  thermocouple  data  to  predictions  showed  some  of  the  dif¬ 
ficulties  in  this  type  of  modeling.  There  are  uncertainties  in  the  meterological 
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Julian  Day 

Figure  12.  Tree  Model  Predictions  for  and  Thermocouple  Data  From  “Joyce 
Kilmer”  for  8  and  9  September  1990.  Data  are  from  the  northern  side  about 
1”  from  the  tmnk  surface 


Julian  Day 

Figure  13.  Tree  Model  Predictions  for  and  Thermocouple  Data  From  “Joyce 
Kilmer”  for  8  and  9  September  1990.  Data  are  from  the  southern  side,  about 
1”  from  the  trunk  surface 
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data,  tree  material  properties,  and  internal  material  distribution.  Moisture  transport 
effects  on  the  model  thermal  response  was  also  not  considered.  Reducing  the  con¬ 
trollable  uncertainties  in  meteorological  data,  refining  the  forced  convective  term, 
and  obtaining  a  better  representation  of  the  material  properties  will  overcome  some 
of  the  model  limitations. 

4.2.2  TVunk  Temperatures  for  Hunter-Liggett  TYee 

A  representative  cross-section  of  the  main  trunk  of  a  mountain  oak  located  at 
Hunter-Liggett,  California  was  modeled.  This  tree  was  selected  as  a  representative 
tree  at  the  site  modeled  for  the  second  year  demonstration  in  the  BTI/SWOE  pro¬ 
gram.  Meteorological  data  for  18  -  20  September  1989  from  this  site  were  used  to 
estimate  the  solar  and  downward  infrared  fluxes.  Figure  16  gives  the  temperatures 
as  a  function  of  time  for  surface  elements  corresponding  to  the  cardinal  compass 
points. 


Figure  16.  Predictions  of  the  Surface  Temperatures  as  a  Function  of  Time  for 
Elements  Near  the  Cardinal  Points.  Calculations  were  made  for  meteorological 
conditions  for  19  -  20  September  1989  at  Hunter-Liggett,  California 

The  impact  of  the  fog  during  the  morning  of  20  September  1989  (see  Table  3) 
on  the  model  thermal  response  can  be  seen  in  a  comparision  of  the  temperature 
predictions  between  the  mornings  of  19  and  20  September  1989. 

On  20  September  1989,  the  temperature  at  the  cardinal  points  tended  to  rise 
at  about  the  same  pace  until  the  fog  burned  off.  This  is  a  result  of  the  way  the 
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model  treats  infrared  and  diffuse  solar  energy.  These  terms  are  implemented  as 
radially  symmetric  to  the  modeled  cylindrical  geometric  shapes.  In  addition,  the 
magnitude  of  these  terms  is  higher  for  the  fog  conditions  of  20  September  1989 
than  for  the  clear  conditions  of  the  previous  morning.  The  dip  in  the  north  and 
west  temperatures  starting  at  about  8  AM  is  a  combination  of  the  fog  cover  burning 
off  and  an  increase  in  the  wind  velocity. 

The  effects  of  clear  conditions  on  19  September  1989  and  solar  position  on 
the  model  temperature  predictions  are  shown  at  the  east  and  west  location.  The 
calculated  east  temperature  rises  in  the  morning  until  that  location  becomes  shaded. 
The  west  location  does  not  receive  direct  solar  energy  until  the  mid-aftemoon.  Its 
temperature  rise  in  the  morning  is  due  to  diffuse  solar  energy,  conduction,  and 
convection  due  to  wind. 

To  demonstrate  the  effect  of  solar  position  on  the  temperature  distribution, 
contour  plots  of  TREETHERM’s  temperature  predictions  through  a  cross  section 
were  made  for  9  PM  and  3  PM  for  20  September  1989.  These  are  shown  in 
Figures  17  and  18.  The  highest  thermal  gradients  are  in  the  location  corresponding 
to  the  solar  vector.  Heating  (or  cooling)  in  the  shaded  section  is  accomplished  by 
the  surface  interaction  with  infrared  radiation,  wind,  and  the  diffuse  component  of 
the  solar  energy. 
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TEMPERATURE  DISTRIBUTION  9  AM.  9/20/89 

0.29  -0.23  -0.17  -0.11  -O.OS  0.01  0.07  0.12  0.18  ( 
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Figure  17.  Model  Temperature  Cross  Sections  for  9  AM,  20  September  1989. 
North  is  to  the  right  and  east  is  at  the  bottom  of  the  figure 

TEMPERATURE  DISTRIBUTION  3  PM,  9/20/89 

-0JI9  -0.23  -0.17  -0.11  -0.05  0.01  0.07  0.12  0.18  OM 


-0.29  -0J3  -0.17  -0.11  -0.05  0.01  0.07  0.12  0.18  0J24 

Figure  18.  Model  Temperature  Cross  Sections  for  3  PM,  20  September  1989. 
North  is  to  the  right  and  east  is  at  the  bottom  of  the  figure 


5  SUMMARY  AND  RECOMMENDATIONS  FOR  FUTURE  WORK 

5.1  Summary 

A  first  generation  3-D  tree  thermal  model,  TREETHERM,  has  been  developed 
for  the  Balanced  Technology  Initiative  on  Smart  Weapons  Operability  Enhance¬ 
ment  Program  for  use  in  scene  simulation  studies.  The  model  can  also  be  used  in 
remote  sensing  and  climate  studies.  TREETHERM  describes  a  single  tree  in  terms 
of  primitive  shapes  and  includes  the  variation  in  material  properties.  The  tree  repro¬ 
duces  the  thermal  response  phenomenology  as  studied  during  a  field  lest  program 
conducted  in  conjunction  with  NASA’s  Forest  Ecosystem  Dynamics  program. 

TREETHERM  predicts  the  overall  thermal  trends  in  an  individual  tree  reason¬ 
ably  well.  The  model  demonstrates  the  thermal  response  of  woody  materials  and 
leaves  to  changes  in  solar  loading,  wind  speed  and  direction,  and  evapotranspira- 
tion.  However,  additional  development  efforts  are  planned  to  refine  and  improve 
the  model.  Specific  efforts  will  include  the  incorporation  of  a  formulation  for  inter¬ 
nal  moisture  flow,  enhancements  to  the  leaf  energy  budget  model,  and  interactions 
with  other  objects  such  as  adjacent  trees  or  other  objects. 

5.2  Recommendations  for  Future  Work 

5.2.1  Addition  of  Woody  Material  Moisture  Transport 

A  full  moisture  treatment  within  the  tree  was  not  included  in  this  version  of 
TREETHERM  due  to  the  complexity  of  describing  the  moisture  content  of  woody 
material  and  its  relationships  to  environmental  conditions.  Instead,  the  moisture 
content  in  the  woody  material  was  approximated  and  the  material  properties  of  dry 
wood  adjusted  and  then  kept  constant  during  the  calculations.  While  acceptable 
for  a  first  order  development  effort,  this  is  not  acceptable  for  a  model  intended  to 
fully  represent  a  tree’s  response  to  changes  in  environmental  conditions. 

The  moisture  content  of  trees  can  change  during  the  course  of  a  day  as  a  result 
of  evapotranspiration  processes  and  in  result  to  changes  in  general  weather  condi¬ 
tions,  such  as  the  presence  of  precipitation.  A  first  order  treatment  could  include 
a  temporal  variation  in  material  thermal  properties.  A  more  detailed  treatment, 
however,  would  require  the  development  of  a  mass  balance  model  for  the  tree. 

5.2.2  Enhancements  to  the  Leaf  Energy  Budget  Model 

The  leaf  energy  budget  model  utilized  a  number  of  simplifications  that  should  be 
readdressed  in  a  second  generation  model.  First,  the  leaf  energy  budget  calculations 
were  performed  for  leaves  assumed  to  be  fully  exposed  to  the  sun.  In  an  actual  tree, 
the  amount  of  sunlight  reaching  inner  leaves  will  be  attenuated  by  the  presence  of 
intervening  leaf  clusters.  The  full  tree  model  has  the  capability  to  keep  track  of 
individual  leaf  clusters  and  their  ability  to  shade  other  portions  of  the  tree.  This 


30 


feature  should  be  coupled  to  the  leaf  energy  budget  model  to  correctly  account  for 
the  solar  loading  on  individual  leaf  clusters. 

Second,  the  convective  exchange  term  needs  to  be  revised  to  account  for  dif¬ 
ferent  shaped  leaves.  The  first  generation  model  assumed  circular  plates  for  the 
leaves.  Other  shapes  should  be  included  to  more  accurately  represent  the  variation 
in  leaf  morphology. 

Finally,  the  evapotranspiration  term  needs  to  be  reexamined.  In  particular,  the 
temporal  and  environmental  variability  in  the  rate  of  water  vapor  diffusion  needs 
to  included. 

5.2.3  TYee-Environmental  Interactions 

This  study  focused  on  the  energy  budget  aspects  of  tree-environmental  inter¬ 
actions.  There  are  other  aspects  of  this  interaction  that  deserve  treatment  in  future 
versions  of  the  model.  One  is  the  effect  of  the  tree  on  wind  fields  near  the  surface. 
In  addition  to  providing  a  source  of  solar  shading  to  the  ground,  trees  can  provide 
a  form  of  wind  shading.  That  is,  the  presence  of  trees  (or  any  3-D  object)  will 
moderate  the  wind  field  near  the  surface  and  this  modification  could  potentially 
have  an  impact  on  the  convective  energy  exchange  at  and  near  the  surface. 

A  second  interaction  deserving  further  study  is  the  thermal  interaction  of  the  tree 
with  the  ground  and  vice  versa.  The  upwelling  ground  infrared  flux  is  a  significant 
term  in  the  tree  and  leaf  energy  budgets.  In  the  present  version,  it  was  not  explicitly 
included  but,  instead,  the  infrared  sources  were  treated  as  coming  from  a  spherically 
symmetric  source.  This  was  done  because  the  inclusion  of  upwelling  infrared 
sources  requires  that  the  tree  be  coupled  to  a  full,  three  dimensional  model  of  the 
surrounding  environment  including  a  treatment  that  accounts  for  shading  effects 
from  the  tree  and  any  other  objects. 
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Appendix  A 
Tree  Shading  Model 

The  tree  shading  model  uses  ray  casting  technology  implemented  for  rendering 
images  or  signatures  of  various  objects  as  observed  by  optical  sensor  systems. 
The  constmctive  solid  geometry  model  in  SENSORSIM  has  been  considerably 
simplified  for  the  current  application,  so  that  trees  are  represented  by  a  series  of 
coaxial  pairs  of  cylinders,  the  inner  cylinder  of  each  pair  representing  a  volume 
of  branch  (or  trunk)  and  the  outer  cylinder  representing  a  volume  filled  to  some 
density  with  leaves.  The  model  determines  the  effective  illuminated  area  for  each 
section  of  the  branch  attenuated  for  partial  shading  by  leaves.  The  model  adaptively 
determines  the  extent  (length  and  width)  and  density  of  the  rays  according  to  the 
location  and  size  of  each  branch  and  leaf  volume.  The  remainder  of  this  section 
describes  the  components  of  this  model  in  greater  detail. 

A-1.  Ray  Casting  Process 

The  ray  casting  process  consists  of  generating  a  grid  of  rays  corresponding  to 
various  lines  of  solar  illumination  (from  the  sun  to  the  tree)  and  determining  what 
part  of  the  tree  is  illuminated  by  each  ray.  The  tree  shading  model  adaptively 
determines  the  extent  of  the  ray  grid  and  the  density  of  rays  in  each  region  of  the 
grid  according  to  the  location  and  size  of  the  cylinders  representing  the  tree.  Since 
the  number  of  rays  which  must  be  cast  is  roughly  proportional  to  the  number  of 
cylinders  used  to  model  the  tree  and  each  ray  must  be  intersected  with  each  cylinder, 
the  execution  time  for  the  shading  model  will  be  approximately  proportional  to  the 
square  of  the  number  of  cylinders  used  to  model  the  tree. 

A  ray  along  the  line  of  solar  illumination  can  be  represented  by  a  point  on  the 
line  (called  the  “ray  origin”)  and  a  direction  vector  parallel  to  the  line  (called  the 
“ray  direction”).*  In  vector  form,  any  ray  is  given  by  the  locus  of  points 

*  Mathematically,  the  terms  “ray”  and  “line”  are  not  synonymous;  a  “line”  ex¬ 
tends  to  infinity  in  both  directions,  while  a  “ray”  begins  at  a  defined  point  and 
extends  to  infinity  in  only  one  direction.  Mathematicians  additionally  apply  the 
term  “half  line”  to  a  ray  with  the  end  point  removed.  All  three  of  these  elements 
can  be  represented  by  Eqn  (A-1),  with  appropriate  restrictions  on  the  value  of  the 
ray  coordinate  /t-  G  3?:  0  <  /.:  <  oo  for  a  half  line,  0  <  A’  <  oo  for  a  ray,  and 
-oo  <  A  <  oo  for  a  line.  In  practice,  rays  originating  at  large  distances  from  the 
target,  and  implementation  considerations  of  processing  precision  dictate  that  the 
ray  origin  should  be  offset  to  a  convenient  location  in  the  proximity  of  the  target. 
As  a  consequence  of  this  offset,  negative  values  of  k  will  ordinarily  occur  and  do 
in  fact  correspond  to  valid  intersection  points;  thus,  the  terms  “ray”  and  “line”  will 


34 


(A-1) 


X{k)  =  Xo  +  A:  AX 

where  X(A:)  is  a  point  on  the  line,  Xq  is  the  ray  origin,  AX  is  the  ray  direction 
vector,  and  k  is  the  ray  coordinate  of  point  X(A:),  The  vectors  in  Eqn  (A-1)  can 


be  separated  into  their  respective  components 

X{k)  =  [x{k)  y{k)  ^(A’)]  (A-2a) 

Xo  ~  [  Vo  ]  (A-2b) 

AX=[Aa:  Ay  Az]  (A-2c) 

to  obtain  the  equivalent  scalar  parametric  equations  for  the  line 

x{k)  =  Xq-^  kAx  (A-3a) 

y{k)  =  yo  +  kAy  (A-3b) 

z{k)  =  2o  +  kAz.  (A-3c) 


The  ray  coordinate  (k)  uniquely  identifies  each  point  along  a  line  of  illumina¬ 
tion,  and  in  fact  is  used  in  the  simulation  to  sort  the  intersected  segments  into  the 
correct  sequence. 

A-2.  Ray  TVansformations 

Before  a  ray  can  be  intersected  with  a  cylinder,  it  must  be  transformed  from 
the  global  (scene)  coordinate  system  into  the  coordinate  system  of  the  cylinder. 
This  transformation  can  be  represented  by  a  square  matrix,  containing  the  magnifi¬ 
cation  and  rotation  information,  which  is  multiplied  by  the  coordinate  vector,  and 
a  translation  vector  which  is  added  to  the  result  of  the  multiplication.  A  point  (or 
any  other  vector  representing  absolute  coordinates)  is  thus  transformed  by* 

Y  =  XA  -I-  B  (A-4) 

where  Y  is  the  transformed  point,  X  is  the  untransformed  point,  A  is  the  magni- 

be  used  interchangeably  in  this  discussion  as  convenience  dictates. 

*  The  choice  of  notation  in  Eqn  (A-4),  coupled  with  the  previous  use  of  X 
(and  the  associated  Xq  and  AX)  for  the  line  of  solar  illumination  in  Eqn  (A-1) 
assume  the  forward  direction  of  all  coordinate  transformations  to  be  from  scene 
coordinates  to  cylinder  coordinates.  Since  both  forward  and  reverse  transformations 
are  available  in  the  simulation,  either  this  or  the  reverse  convention  could  be 
selected  arbitrarily;  for  consistency  of  notation,  this  convention  will  be  retained. 
If  the  reverse  convention  had  been  chosen  (that  is,  if  the  forward  direction  of  a 
transformation  mapped  from  primitive  shape  coordinates  to  scene  coordinates),  the 
effect  would  have  been  to  interchange  X  and  Y  in  Eqns  (A-4),  (A-5),  and  (A-6), 
and  to  interchange  the  forward  transformation  matrices  A  and  B  with  the  reverse 
transformation  matrices  A'  and  B'  respectively  in  all  subsequent  equations. 
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Table  A-1.  Primitive  Transformations 


Transformation  A  B 


ri 

0 

0] 

ri 

0 

01 

TVanslation 

0 

1 

0 

[  ~iy  ~^z  ]  0 

1 

0 

[<X  ty  <x] 

.0 

0 

1. 

.0 

0 

1. 

■1/sx  0  0  I  px  0  O' 

Scaling  0  l/sj,  0  [0  0  0]  0  0  [0  0  0] 

.  0  0  1/sJ  Lo  0  Ax. 

■  1  0  0  1  r  ^  °  ' 

Rotation  -  X  Axis  0  cosip  —  sin^'  [0  0  0]  0  cosip  sin  0  [0  0  0] 

.0  sinip  cos  Ip  J  Lo  —  sin  ^  cc^ip, 

cos  V'  0  sin  ^1  r  cos  ip  i)  —  sin  V’ " 

Rotation  -  Y  Axis  0  10  [0  0  0]  0  10  [0  0  0] 

.  — sinV'  0  cosV’J  LsinV*  0  cosip 

'cosip  —sinip  O']  TcosV’  sin^l>  0" 

Rotation  -  Z  Axis  sin^'  cos  ip  0  [0  0  0]  sin^  cosip  0  [0  0  0] 

.  0  0  ij  L  0  0  ij 

fication  and  rotation  matrix,  and  B  is  the  translation  vector.  The  inverse  of  this 
transformation  is  obtained  by  solving  Eqn  (A-4)  for  X 


X  =  YA-l-BA-l  (A-5) 

The  inverse  transformation  can  thus  be  represented  by 

X  =  YA'  +  B'  (A-6) 

where 

A'  =  A-1  (A-7a) 

B'  =  -BA-1.  (A-7b) 


The  values  of  A,  A',  B,  and  B'  for  the  primitive  transformations  are  shown  in 
Table  A-1. 

A  mle  for  combining  two  transformations  may  easily  be  derived  by  letting 
subscripts  1  and  2  be  applied  to  the  transformation  matrices  to  designate  the  two 
transformations  to  be  combined,  and  the  unsubscripted  matrices  designate  the  com¬ 
bined  transformation.  The  forward  and  reverse  composite  transformations  are  given 
by 

Y  =  (XAi-hBi)A2-hB2 

X  =  (yA'2  +  bQ  A'l  +  B'l 
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(A-8a) 

(A-8b) 


Performing  the  matrix  algebra  necessary 
(A-6), 

to  obtain  the  form  of  Eqns  (A-4)  and 

A  =  A1A2 

(A-9a) 

A'  =  A'2A' 

(A-9b) 

B  =  B1A2  -f  B2 

(A-9c) 

B'  =  B^A'i4-Bi.' 

(A-9d) 

Tliis  derivation  makes  no  assumptions  about  the  transformations  being  combined, 
so  these  mles  may  be  used  recursively  to  compute  arbitrary  combinations  of  trans¬ 
formations.*  In  the  tree  shading  model,  three  transformations  are  applied  for  each 
cylinder,  consisting  of  (1)  translating  the  coordinate  origin  to  the  center  of  the 
cylinder,  (2)  rotating  the  coordinate  system  to  coincide  with  the  axes  of  the  cylin¬ 
der  (using  a  precomputed  rotation  matrix  composed  of  several  primitive  rotations), 
and  (3)  scaling  to  ^e  cylinder  to  the  proper  dimensions  (the  horizontal  (x  and 
y)  components  are  scaled  by  half  the  radius  of  the  cylinder,  and  the  vertical  (z) 
component  is  scaled  by  half  of  the  height). 

A-3.  Ray-Cylinder  Intersections 

The  untransformed  cylinder  has  a  radius  of  one  unit  and  a  height  of  two  units, 
and  is  centered  at  the  origin  of  its  local  coordinate  system.  The  axis  of  the  cylinder 
is  the  vertical  (z)  axis.  This  cylinder  has  a  volume  of  Stt  cubic  meters,  and  is 
characterized  by  the  inside-outside  function 

f{x,  y,  z)  =  max  -1-  -  l)  ,  2  -  1,  -1  -  2)  .  (A-10) 

The  boundary  of  the  cylinder  is  the  locus  of  points  where  the  inside-outside  function 
evaluates  to  zero. 

The  incident  ray  can  be  expressed  by  separating  Eqn  (14)  into  its  respective 
components 

x{k)  =  Xo  +  kAx  (A-1  la) 

y{k)  =  yo  +  kAy  (A-1  lb) 

*  The  ability  to  arbitrarily  combine  transformations  by  this  mle  is  actually  de¬ 
pendent  on  the  associative  property  of  the  operation  of  combining  transformations 
defined  by  Eqns  {A  —  9a)  through  {A  —  9d).  This  property  can  be  proven  by 
using  the  associative  and  distributive  laws  of  matrix  multiplication  and  addition 
to  show  that  A  =  (AiA2)A3  =  Ai(A2A3)  and  B  =  (B1A2  -f  B2)A3  4-  B3  = 
Bi(A2A3)  -f  (B2A3  +  B3).  It  is  not  necessary  to  prove  the  corresponding  rela¬ 
tionships  for  the  inverse  transformations  because  the  inverse  transformations  are 
themselves  valid  transformations  and  thus  governed  by  the  same  theorem. 


37 


(A- 11c) 


z{k)  =  2o  +  kAz 


where 


Yq  —  [a^'o  Vo  ^o] 
AY=[Ax  Ay  Az] 


(A- 12a) 
(A- 12b) 


The  points  where  the  ray  intersects  the  surface  of  the  cylinder  may  be  obtained  by 
first  considering  the  infinite  cylinder  represented  by  the  first  argument  to  the  maxi¬ 
mum  function  in  Eqn  {A  —  10).  Setting  this  expression  to  be  zero  and  substituting 
Eqns  (.4  —  11a)  and  {A  —  116)  yields  an  equation  for  the  intersection  of  the  ray 
with  the  infinite  cylindrical  surface  of  the  form 


ak"^  -\-bk  A  c  =  0 


(A-13) 


where 


a  =  Ax^  -1-  Ay^  (A-14a) 

6  =  2(xoAx  +  yoAy)  (A-14b) 

c  =  Xo-fyo-l.  (A-14c) 

If  a  is  not  zero,  the  roots  of  Eqn  (13)  are  given  by  the  quadratic  formula 


-(6/21 


{b/2r—ac 


(A- 15) 


The  intersection  points  given  by  Eqn  {A  —  15)  are  complex  or  the  same  whenever 
the  ray  misses  (that  is,  does  not  intersect)  the  infinite  cylinder.  This  condition  is 
identified  by  the  discriminant  of  the  quadratic  (the  quantity  under  the  square  root  in 
Eqn  {A  -  15)).  If  the  discriminant  is  negative  or  zero,  the  ray  misses  the  cylinder 
and  no  further  computation  is  performed;  otherwise,  the  ray  intersects  the  infinite 
cylinder  in  exactly  one  segment.  Since  a  cannot  be  negative  (it  is  the  sum  of  the 
squares  of  two  real  values)  and  the  special  case  where  a  is  zero  requires  separate 
treatment,  the  ray  coordinate  corresponding  to  the  smaller  range  (the  point  where 
the  ray  enters  the  infinite  cylinder)  is  obtained  by  choosing  the  negative  sign  for 
the  square  root,  and  the  ray  coordinate  corresponding  to  the  greater  range  (where 
the  ray  exits  the  infinite  cylinder)  is  obtained  by  choosing  the  positive  square  root. 

If  the  ray  is  found  to  intersect  the  infinite  cylinder,  it  is  a  simple  matter  to 
verify  whether  it  also  intersects  the  finite  cylinder.  To  perform  this  verification, 
the  vertical  (2)  component  of  the  infinite  cylinder  entrance  and  exit  points  are 
computed  by  substituting  the  respective  roots  obtained  from  Eqn  {A  —  15)  into 
Eqn  (.4  -  lie).  If  both  points  are  at  or  above  the  top  of  the  finite  cylinder  (that  is, 
2  >  1  for  both  roots),  the  ray  misses  (that  is,  does  not  intersect)  the  finite  cylinder. 
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The  ray  likewise  misses  the  finite  cylinder  if  both  points  are  at  or  below  the  bottom 
of  the  finite  cylinder  (that  is,  ^  <  - 1  for  both  roots).  In  either  of  these  cases, 
no  further  computation  is  necessary.  In  all  other  case,  the  ray  intersects  the  finite 
cylinder  in  some  subset  of  its  intersection  with  the  infinite  cylinder. 

It  exactly  one  endpoint  of  the  segment  where  the  ray  intersects  the  infinite 
cylinder  is  above  the  top  of  the  finite  cylinder,  the  corresponding  endpoint  of  the 
intersection  segment  for  the  finite  cylinder  must  be  on  the  top  of  the  finite  cylinder. 
The  vertical  (c)  coordinate  of  this  intersection  point  is  that  of  the  top  of  the  cylinder 
(::  =  !,  obtained  by  setting  the  second  argument  of  the  maximum  function  in  Eqn 
{A  -  10)  to  zero);  combining  this  result  with  Eqn  {A  -  11c)  and  solving  for  the 
ray  coordinate 


k  = 


(A- 16) 


Likewise,  if  exactly  one  endpoint  of  the  segment  where  the  ray  intersects  the 
infinite  cylinder  at  a  point  which  is  below  the  bottom  of  the  finite  cylinder,  the 
corresponding  endpoint  of  the  segment  where  the  ray  intersects  the  finite  cylinder 
is  given  by 


k  = 


(A-17) 


It  is  not  necessary  to  consider  the  special  case  of  zero  values  of  in  Eqns  (.4-16) 
and  {A  —  17),  since  this  condition  would  correspond  precisely  to  a  horizontal  ray 
which  could  neither  enter  nor  exit  through  the  top  or  bottom  of  the  cylinder. 

The  case  of  a  vertical  ray  (a  =  0)  requires  special  treatment,  as  Eqn  (A  -  13) 
becomes  irrelevant.  (By  Eqn  (A  -  14tt),  a  is  zero  if  and  only  if  both  Ax  and  Ay 
are  zero,  in  which  case  b  is  also  zero  by  Eqn  (A  -  146).)  In  this  case,  the  ray 
intersects  the  infinite  cylinder  if  and  only  if  c  is  negative;  thus,  if  c  is  not  negative, 
no  further  computation  is  necessary.  If  c  is  negative,  the  intersection  points  must 
be  the  top  and  bottom  of  the  cylinder,  and  the  ray  coordinates  can  be  computed 
directly  from  Eqns  (A  -  16)  and  (A  -  17).  A  check  of  the  sign  of  Az  is  sufficient 
to  tell  whether  the  ray  enters  the  cylinder  at  the  top  and  exits  at  the  bottom  (A^ 
negative)  or  vica  versa  (Az  positive);  Az  cannot  be  zero  since  the  ray  direction 
vector  cannot  be  a  zero  vector. 


The  surface  orientation  vector  (normal  to  the  cylinder)  in  primitive  shape  co¬ 
ordinates  is  obtained  by  combining  Eqn  (A  —  10)  with  the  relationship 


Ny  =  V/ 


(A- 18) 


For  a  ray  intersecting  the  cylindrical  part  of  the  surface,  the  first  argument  to 
the  maximum  function  in  Eqn  (A  —  10)  prevails  and  the  normal  vector  is  given  by 

Ny  =  [  a;  y  0  ] 
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(A- 19) 


wherea:  and  y  are  computed,  respectively,  from  Eqns  {A  — Ha)  and  (A— 116).  For 
segment  endpoints  on  the  top  of  the  cylinder,  the  second  argument  to  the  maximum 
function  prevails  and  the  normal  vector  is 

Ny  =-  [  0  C  1  ]  (A-20) 

Likewise,  on  the  bottom  of  the  cylinder,  tlie  third  argument  to  the  maximum 
function  prevails  and  the  normal  vector  is 

Ny  =  [0  0  -1]  (A-21) 

Since  the  cylinder  is  defined  to  be  a  mathematical  open  set,  the  vertical  coordinate 
of  the  intersection  point  can  be  used  to  select  the  correct  equation  for  the  normal 
vector.  The  normal  vector  is  used  to  determine  the  sector  of  the  branch  which  is 
intersected  by  the  ray. 

A-4.  Effective  Illuminated  Area 

The  effective  illuminated  area  is  defined  to  be  the  area  of  a  planar  patch  oriented 
normal  to  the  direction  of  solar  radiation  that  would  receive  the  same  incident  solar 
radiation  as  a  sector  of  a  branch  (or  trunk).  To  compute  the  effective  illuminated 
area  for  each  sector  of  each  branch,  an  effective  area  for  each  ray  is  computed  by 
multiplying  the  area  corresponding  to  each  ray  (the  square  of  the  ray  separation  in 
the  applicable  region  of  the  ray  grid)  by  any  attenuation  due  to  partial  shadowing 
by  leaf  volumes  intersected  before  the  first  intersected  branch  cylinder,  which  is 
assigned  to  the  appropriate  sector  of  the  first  branch  cylinder  intersected  by  the 
ray.  The  effective  area  for  rays  incident  upon  the  ends  of  cylinders  are  neglected. 
For  rays  intersecting  the  cylindrical  surface,  the  intersected  sector  is  determined 
from  the  orientation  of  the  surface  normal  vector  in  the  local  coordinate  system. 
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